Milestone 6: The Sliding Lid

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as anim
import IPython.display as display
import src.simulation as sim
import src.boundary as bdry
import src.visualization as visz
import os
import sys
from tqdm import trange, tqdm

Parameters

Combination of 12 different parameter sets ($L$, $\nu$, $U_w$) and the resulting Reynolds number $Re$:

$L$	$\nu$	$\omega$	$U_w$	$Re$
300	0.03	~1.69	0.025	250
			0.05	500
			0.075	750
			0.1	1000
	0.09	~1.30	0.075	250
			0.15	500
			0.225	750
			0.3	1000
	0.15	~1.05	0.125	250
			0.25	500
			0.375	750
			0.5	1000

L = 300
n_steps = 100000
params = [(Re, nu, Re*nu/L)
          for nu in [0.03, 0.09, 0.15]
          for Re in [250, 500, 750, 1000]]

# animation: plot only every nth step and increase step size over time
anim_stages = visz.PlotStages([(0, 10), (1000, 50), (5000, 100), (20000, 500), (50000, 1000)])
anim_steps = anim_stages.list(n_steps)

Initialize

lattice = [
    sim.LatticeBoltzmann(L, L, viscosity=viscosity, init_density=np.ones((L, L)),
                         boundaries=[bdry.MovingWallBoundaryCondition("t", [0, wall_velocity]),
                                     bdry.RigidWallBoundaryCondition("blr")],
                         plot=False, animate=False)
    for _, viscosity, wall_velocity in params
]

Simulation and Measurements

for l in range(len(lattice)):
    velocity_field = np.empty((len(anim_steps), 2, L, L))
    j = 0
    for i in trange(n_steps, file=sys.stdout, desc=f"lattice {l+1:2}/{len(lattice)}"):
        if i in anim_steps:
            velocity_field[j] = lattice[l].velocity
            j += 1
        lattice[l].step()
    velocity_field[-1] = lattice[l].velocity
    # save result
    os.makedirs("out/m6", exist_ok=True)
    np.save(f"out/m6/lattice_{l}.npy", velocity_field)

lattice  1/12: 100%|██████████| 100000/100000 [33:45<00:00, 49.38it/s]
lattice  2/12: 100%|██████████| 100000/100000 [33:16<00:00, 50.10it/s]
lattice  3/12: 100%|██████████| 100000/100000 [34:24<00:00, 48.44it/s]
lattice  4/12: 100%|██████████| 100000/100000 [33:15<00:00, 50.11it/s]
lattice  5/12: 100%|██████████| 100000/100000 [33:08<00:00, 50.30it/s]
lattice  6/12: 100%|██████████| 100000/100000 [33:18<00:00, 50.04it/s]
lattice  7/12: 100%|██████████| 100000/100000 [33:11<00:00, 50.21it/s]
lattice  8/12: 100%|██████████| 100000/100000 [33:24<00:00, 49.89it/s]
lattice  9/12: 100%|██████████| 100000/100000 [33:05<00:00, 50.38it/s]
lattice 10/12: 100%|██████████| 100000/100000 [33:07<00:00, 50.31it/s]
lattice 11/12: 100%|██████████| 100000/100000 [32:23<00:00, 51.45it/s]
lattice 12/12: 100%|██████████| 100000/100000 [33:14<00:00, 50.15it/s]

Plot Results

Steady-state streamplots for different parameters

n_cols = 4
n_rows = len(lattice) // n_cols
fig, ax = plt.subplots(n_rows, n_cols, sharex=True, sharey=True, figsize=(19, 16*(n_rows/n_cols)), gridspec_kw={"width_ratios": [1, 1, 1, 1.11]})
axf = ax.flatten()
axf[0].set_xlim(0, L)
axf[0].set_ylim(0, L)
axf[0].invert_yaxis()
x = np.arange(L)
y = np.arange(L)
for l in range(len(lattice)):
    Re, viscosity, wall_velocity = params[l]
    axf[l].set_title(f"$Re={Re}$, $\\nu={viscosity}$, $\omega \\approx {lattice[l].omega:.1f}$, $U_w={wall_velocity}$")
    if l % n_cols == 0:
        axf[l].set_ylabel("$y$ dimension")
    if l >= n_cols*(n_rows-1):
        axf[l].set_xlabel("$x$ dimension")
    velocity_field = np.load(f"out/m6/lattice_{l}.npy")
    u = velocity_field[-1, 1]
    v = velocity_field[-1, 0]
    norm = plt.Normalize(0, wall_velocity)
    axf[l].streamplot(x, y, u, v, color=np.sqrt(u**2 + v**2), norm=norm, density=1.2)
    if l % n_cols == 3:
        cbar = fig.colorbar(plt.cm.ScalarMappable(norm, plt.cm.viridis), ax=axf[l], fraction=0.05,
                            ticks=[0, wall_velocity], label="velocity magnitude $|u|$")
        cbar.ax.set_yticklabels(["0", "$U_w$"])
fig.tight_layout()
plt.show()

Animation

%matplotlib ipympl

l = 7
Re, viscosity, wall_velocity = params[l]
velocity_field = np.load(f"out/m6/lattice_{l}.npy")
velocity_magnitude = np.sqrt(velocity_field[:, 1]**2 + velocity_field[:, 0]**2)

fig, ax = plt.subplots(figsize=(6.8, 6.0), dpi=180)
fig.suptitle(f"Sliding Lid Experiment", fontweight="bold")
cmap = plt.cm.viridis
fig.colorbar(plt.cm.ScalarMappable(plt.Normalize(0, wall_velocity), cmap), ax=ax, label="velocity magnitude $|u|$")
x = y = np.arange(L)


def update_streamplot(frame):
    i, t = frame
    ax.clear()
    ax.set_title(f"$Re={Re}$, $\\nu={viscosity}$, $\omega \\approx {lattice[l].omega:.1f}$, " +
                 f"$U_w={wall_velocity}$, $t={t}$ (x{anim_stages.get_step_size(t)})", loc="left")
    ax.set_xlabel(f"$x$ dimension")
    ax.set_ylabel(f"$y$ dimension")
    ax.set_xlim(0, L)
    ax.set_ylim(0, L)
    ax.invert_yaxis()
    fig.tight_layout()
    u, v = velocity_field[i, 1], velocity_field[i, 0]
    c = velocity_magnitude[i] if velocity_magnitude[i].sum() > 0 else "black"
    return ax.streamplot(x, y, u, v, density=2.0, color=c, cmap=cmap)


animation = anim.FuncAnimation(fig, update_streamplot, frames=list(enumerate(anim_steps)))
with tqdm(anim_steps, file=sys.stdout) as pbar:
    animation.save("out/m6/sliding_lid.webm", fps=25,
                   writer="ffmpeg", codec="libvpx-vp9", extra_args=["-an", "-crf", "42"],
                   progress_callback=lambda i, n: pbar.update())
plt.close()

100%|██████████| 441/441 [06:15<00:00,  1.17it/s]

display.Video("out/m6/sliding_lid.webm", embed=True, height=540)